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I. INTRODUCTION 

Much attention has been focused on Bose-Einstein con- 
densation (BEC) realized in atomic gases |IJ Q pi Q. 
The atomic species investigated include the isotopes ° 7 Rb 
and 85 Rb and also 23 Na, 7 Li, H, 4 He, and 41 K (see 
Refs. § |, 1 1, ||, 0, §). These systems provide one 
with unique opportunities to investigate novel states of 
superfluid matter waves. When an atomic gas is cooled 
in an applied external magnetic field, the condensate is 
described by a scalar order parameter. The magneto- 
optic trapping method is widely used to study BEC. Re- 
cently it has become possible to trap an atomic gas using 
purely optical methods, thus the resulting BEC retains 
its original atomic hyperfine state. Specifically, 23 Na j|| 
and 87 Rb both with F = 1 are successfully Bose- 
condensed via optical trapping. This system, dubbed 
spinor BEC, is now described by multicomponent order 
parameter. 

Ohmi and Machida JO], and Ho jl2| have indepen- 
dently introduced the basic Hamiltonian for describing 
this spinor BEC by extending the Bogoliubov formal- 
ism to study the fundamental properties of this interest- 
ing multicomponent BEC, pointing out the richness of 
the topological defect configurations such as the Z-vector 
textures and domain-wall structures. These are analo- 
gous to those found in the A and B phases of superfluid 
3 He H, 11- An advantage of the dilute BEC sys- 
tems with respect to strongly interacting liquid helium 
(for 4 He, there exists no microscopic theory; for 3 Hc, 
there is a BCS-like microscopic theory, generalized for 
p-wave pairing) is that for the dilute gas one can make 
controlled approximations, treating the interparticle in- 
teractions as small pertubative parameter. Moreover, it 
is possible to directly visualize the condensate in atomic 
gases using optical methods. Spin of condensate is easily 
controlled by external magnetic field. The BEC systems 
are quite versatile also in that the interaction parameter 
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can be adjusted over a large range — even its sign can 
be changed. Finally, the condensate in BEC systems can 
exhibit several spin states, e.g., with the hyperfine spins 
F = I and F = 2. 

Theoretical studies of vortices and other topologi- 
cal defects in spinor BECs was initiated by Ohmi and 
Machida |[f| and Ho jl^j . Systematic investigations 
on vortices were followed by Yip |fl6f who considered 
both axisymmetric and non-axisymmetric vortices, and 
by Isoshima et al. p7| , who only considered axisymmet- 
ric vortices and their excitation spectra. Leonhardt and 
Volovik Jl|, Stoof fH, Merzlin et al. @, Zhou jp and 
Martikainen and Suomincn p2[ examined exotic topolog- 
ical defect structures in spinor BECs. 

Here we continue our studies of spinor BECs g, 
f24| , [25| , |26f ; in particular, those on vortices with 
F = 1 [I7j ]. We investigate the vortex phase diagram in 
the plane spanned by the magnetization (which is given 
a priori when a three-component atomic gas is prepared) 
and the external rotation. We consider both the antifer- 
romagnetic and ferromagnetic cases. The former is real- 
ized in 23 Na while the latter is expected for 87 Rb. The 
nonmagnetic situation in which the spin channel interac- 
tion g s = is studied as the limiting case for \g s \ <C g n 
(g s and g n are the interaction constants for the spin and 
density channels). We restrict our calculations to ax- 
isymmetric vortices with winding numbers less than or 
equal to unity |ljj . 

The organization of this paper is as follows: After giv- 
ing a brief introduction to the Hamiltonian for the system 
and the extended Gross-Pitaevskii equation, we enumer- 
ate the possible vortices allowed by axisymmetry in Sec. 
2. To investigate the global stability of the various vor- 
tex structures, the relative free energies of the different 
vortices are compared as functions of magnetization and 
rotation frequency in Sec. 3. Section 4 presents the exci- 
tation spectra for each vortex by solving the associated 
Bogoliubov equations extended to an order parameter 
with three components in order to investigate whether a 
given vortex state is stable against collective modes. This 
yields a local stability criterion for each vortex type. This 
consideration extends our previous works of nucleation 
criteria |2^, E8l E9L BO, pl|. The final Section presents a 



2 



summary and discussion. 



II. POSSIBLE TYPES OF AXISYMMETRIC 
VORTICES 

We treat the system of a Bose condensate with inter- 
nal degrees of freedom F = 1. Hence the condensate or- 
der parameter is characterized by three components with 
nip = 1,0,-1. External rotation of the system around 
the rotation axis, perpendicular to the disc-shaped two- 
dimensional plane, is denoted by the angular velocity f2 
which has a sense (+ or — ). 

We start with the system Hamiltonian jn], [ljj 
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with C = h /(2m a ). The interaction is character- 
ized by the two kinds of channels; the density chan- 
nel: g n = 4nh 2 (ao + 2a2)/(3m a ), and the spin chan- 
nel: g s — 4Tth 2 (ci2 — cio)/(3m a ). The atomic mass is 
m a . The scattering lengths ao and 0,2 characterize colli- 
sions between atoms with total spin and 2. The sub- 
scripts are a — (x, y, z) and i, j, k, I = (0, ±1). The latter 
correspond to the above three species. The scalar field 
V(r) is the harmonic confining potential. The angular- 
momentum matrices F a are 
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F x = 



F„ 



F, = 



The chemical potentials /ij obey /ix — A«o = Mo — M-i- 
We introduce fi = (j,q and \i! = \i\ — (j,q. The Gross- 
Pitaevskii equation for this system, extended to the three 
components becomes 




E 



{-CV 2 + V(r)-(u + »>j)}6 jk 

+gn^2\<t>i\ 2s jk 



lm<Pl Vm 



lh 
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The total energy of the condensate is given by 
E = I ' d 2 r 
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where 
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Since we treat a cylindrically symmetric disc-shaped 
two-dimensional system, the condensate wavefunction <j>j 
may be decomposed in the form: 



Mr, 6) = &(r) 7j -(0) 

= <pj (r) exp[i(aj + (3jt 
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where the condensate wavefunctions </>j(r) (j — 0,±1) 
are expressed in terms of cylindrical coordinates. The 
phases of the condensate wavefunctions are determined 
such that the energy E s in Eq. (||) is minimized. The 
minimizing of the above spin-dependent term: 
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leads to the condition 
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where the upper (lower) sign is used for the ferromagnetic 
(antiferromagnetic) case g s < (g s > 0), respectively. 
This condition may be rewritten in terms of a and (3 as 
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where n is an integer. We take a±i — ag = in the fol- 
lowing since they have no effect in the discussion below. 
The phases of the three components are now expressed 
in the form 
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using j3 = Pi — flo = /3q — In Fig. |], the various 

lines show the possible combinations of the [3 il giving 
rise to the vortex types which we consider in the follow- 
ing. If we restrict our consideration to the phase coef- 
ficient, or the winding number being less than 2, there 
are nine possible combinations of the Pi. There is one 
vortexfree configuration and the other combinations are 
(1,0,-1), (1,1,1), and (1, i, 0) where the winding numbers 
of the three components (pi, (f>o, and 0_i are denoted in 
this order. The latter vortex configuration is only re- 
alized in a form (l,x,0) where cpo vanishes (where Pi is 
fractional). It should be noted that the restriction for the 
winding number to be less than 2 is reasonable, because 
a vortex with a higher winding number is unstable and 
easily breaks up into vortices with unit winding number. 




FIG. 1: Possible allowed vortex configurations under axisym- 
metry. Each condensate 0; (i = 1,0,-1) has the phase /3i. 
The vortices; (1,0,-1), (l,x,Q) and (1,1,1) are treated in this 
paper. 



III. RELATIVE STABILITY OF VARIOUS 
VORTICES 

We compare the relative energies of the various vor- 
tices enlisted above. There are eight candidates for the 
possible vortex type, allowed by axisymmetry. Three of 
them are (1,0,-1), (1,1,1), and (l,x, 0) where "x" denotes 
that this component vanishes identically as mentioned 
above. The last one is named Alice vortex after the Alice 
string p!l[ . By interchanging the roles of <pi and ef>-i, we 
can obtain two more vortices: (-1,0,1), and (0,x,l). The 
remaining combinations (0,a;,l),(-l,a;,0) and (-1,-1,-1) al- 
ways have higher energies than those of (0,x,l), (l,x,0) 
and (1,1,1), respectively, when the external angular ve- 
locity Q is positive. Thus, in what follows, we consider 
the vortices, (1,0,-1), (1,1,1), (l,x,0), (-1,0,1), and (0,x,l) 
for n > and M/N > 0. 

The mass of a 87 Rb atom is m a = 1.443 x 10~ 25 kg. 
This value is used in the following. We use the scatter- 
ing lengths ao = 5.5 x 10 _9 m and 02 = 5.843 x 10~ 9 m 



for the antiferromagnetic case. The ratio of the in- 
teraction coefficients is g s /g n = +0.02. Another set 
for the ferromagnetic case is ao — 5.5 x 10~ 9 m and 
02 = 5.182 x 10~ 9 m, whose ratio of the interaction co- 
efficients is g s /g n = —0.02. The coefficient of the mag- 



netic term vanishes when ao 
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5.5 x 10- y m. The 



ratio becomes g s /g n — 0. The amplitude of the ratio 
\g s /g n \ — 0.02 is for a Na atom. 

For each vortex type and the g s and //, the chemical 
potential /x is determined such that the linear density of 
the particle number becomes 2 x 10 3 (/im) _1 . We use a 
harmonic potential with v = 200 Hz for radial confine- 
ment. All the energies are scaled by the trap frequency 
hv. The angular velocity Q is normalized by 2iri> below. 

The spatial profile of the condensate wavefunctions in 
each vortex is determined by the above GP Eq. (||), which 
is expressed as 
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When we vary fl and //, the condensate wavefunctions 4>i 
are determined by a single parameter TiSip + fi' because 
the variation of htlPo is canceled by that of the chemical 
potential /1. 

The energy of the system Eq. (Q) is rewritten as 
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where the number of particles in the j'th component is 
Nj = J d 2 r\cf>j\ 2 . The total particle number in the sys- 



tem is N = Nj , and the total angular momentum is 
L = + jP)Nj. The energy E ilmet , the magneti- 

zation M = ^2jjNj, and L are functions of the <f>j and, 
therefore, are functions of ft£ip + p! . The critical angular 
velocity Q c for a given magnetization between two vortex 
types (1) and (2) is simply determined by 



hn c (M) = 



E 



(i) 



-K 
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The range of M which a vortex state can take depends 
on the interaction constant g s and the winding number. 
For example, the (1,1,1) vortex does not have an inter- 
mediate value of M/N in the cases with g s — (non- 
magnetic case) and — 0.02<7„ (ferromagnetic case). This 
behavior is well known for the uniform spinor conden- 
sate [Fy], 0|. When g s — —0.02g n , the magnetization 
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range —1 < M/N < is not allowed for the (l,x,0) 
vortex. This range may be improved with improved nu- 
merical accuracy |^9|, [52|. 

Figure 2 shows the phase diagrams with the lowest 
energy in the plane of M/N and SI for the three cases; 
the antiferromagnetic (g s = 0.02g n ), the nonmagnetic 
(g s = 0), and the ferromagnetic cases (g s = — 0.02g„). 
The critical angular velocity O c (M) in Eq. Jl4| ) deter- 
mines the boundaries of each domain. When the system 
is antiferromagnetic or nonmagnetic, the (0,2,1) vortex 
has the lowest energy around SI = for M/N = 1. This 
is because the (0,2,1) vortex almost reduces to the vor- 
texfree single-component state without the vortex wind- 
ing in the full polarization limit. 

The (1,0,-1) vortex is lowest in energy around SI = 
with M/N = 0, irrespective of the magnetism. The con- 
densate simply becomes a vortexfree system for the ferro- 
magnetic case. For the antiferromagnetic case, the nonro- 
tating (f> component without winding number dominates 
and the 4>\ and 0_i components which have opposite 
winding numbers are small. 

This (1,0,-1) vortex and the Alice vortex compete in 
the region around SI ~ 0.6, M/N ~ 0.8. The curve 
between them reflects the kinetic energy of the (1,0,- 
1) vortex when there is no magnetic term (g s = 0). 
The amplitudes of the condensate wavefunctions 4>i(r) 
and (f>o(r) of the (1,0,-1) vortex for the magnetization 
M/N(> 0) are equal to 4>i(r) and 4>-i(r) of the (1,2,0) 
vortex for the magnetization 2M /N— 1. The (1,0,-1) and 
the (1,2,0) vortex have the angular momenta L = M and 
L = (M + N)/2, respectively. The critical Sl c in Eq. |L4|) 
becomes 



Ml c {M) = 2 



E[ lmcr {M) - E[ lm J(M + AQ/2) 
M-N 



(15) 



where E' inneT (M) is i?i rm cr of the (1,0,-1) vortex for < 
M < N. The behavior of Sl c is still complicated because 
both the numerator and the denominator approach zero 
when M approaches N. 



IV. LOCAL STABILITY 

There is another stability criterion for a vortex. The 
excitation spectrum for a stable vortex must be positive 
definite. In other words, if the lowest excited state in the 
spectrum becomes negative as a function of SI, the given 
vortex becomes locally unstable. 



A. Excitation spectrum 

The excitation spectrum is obtained by solving the 
following Bogoliubov equation extended to the three- 
component BEC case: 
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FIG. 2: Phase diagrams in the Q vs M/N plane, showing the 
vortex type which has the lowest energy E(Q). (a) Antiferro- 
magnetic case (<?,, = 0.02<7„). (b) Nonmagnetic case (g 3 — 0). 
(c) Ferromagnetic case (g a = —0.02g„). The (1,1,1) vortex is 
absent from the phase diagrams (b) and (c) because this vor- 
tex cannot have a stable configuration for — 1 < M/N < 1. 
The (0,2,1) Alice vortex is not stabilized and it is thus ex- 
cluded from the phase diagram (c). 
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FIG. 3: Local stability regions bounded by for the two 
types of vortices, (1,0,-1) and (l,x,0). (a) Antiferromagnetic 
case (g B = 0.02g„). (b) Nonmagnetic case (g s = 0). (c) 
Ferromagnetic case (g s = — 0.02<7„). The solid line and the 
dashed line show the upper critical an d the lower critical 
f2~ for a (1,0,-1) vortex. The dashed-dotted line and the 
dotted line show the upper critical f2+ and the lower critical 
Qc f° r the (l,i,0) Alice vortex. The lower critical f2 is higher 
than the upper critical Q for -0.5 < M/N < 0.5 in the (1,0,- 
1) vortex. The range -1 < M/N < -0.1 for the (l,x,0) 
vortex is inhibited because of the ferromagnetic interaction. 
When the imaginary part of e is finite, its real part is used 
in Eqs. ( p^ ) and (|25|). These regions are indicated with the 
crosses and pluses. 
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where u q (r, i) and v q (r, i) are the cj'th eigenfunctions with 
the spin component i and e q corresponds to the g'th 
eigenvalue. The normalization condition for the u q {r, i) 
and v q (r, i) is 



£ J {K(r,*)l 2 -K(M)l 2 }dr 



= i. 



(19) 



In our disc-shaped system, the wavefunctions have 
the phase u q (r,j) — u q (r, j)e l9i - qe+l3o+ ^" > and v q (r,j) — 
v q (r,j)e l0( -i<>-l 3 «-J( 3 l Then Eqs. (|l6| - |lS|) are written as 
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The excitation energy e q varies as a function of the 
angular velocity f2 as determined by Eqs. ( |l6| ) and (|2j 
The critical values for local stability are defined by 



Q+ = mm(e q /q 9 ) 
qe>0 

Q- = max(e q /q e ). 

qe<0 



(24) 
(25) 



The critical velocity f2+ defined above corresponds to 
the instability of the surface excitations; namely, the en- 
ergy of the excitation modes for qg > become negative 
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(see the paper 1 28 by Isoshima and Machida for the one- 
component BEC). The critical velocity fl~ corresponds 
to the local instability where the core excitation modes 
with qg < become negative upon varying f2. In the fol- 
lowing, we determine f2j, yielding a stability region in 
the f2 vs M/N plane where the all excitation modes are 
positive definite for a given vortex configuration. Namely, 
we evaluate flf for the vortex types (1,0,-1), (1,1,1), and 
(1,2,0) as functions of fl and M/N. 

When there is a mode with qg = and e < 0, the 
system has an instability regardless of external rotation 
and the definition of flf becomes meaningless. There are 
few modes found with e < for qg — 0. The value is at 
— 10 -12 < e < 0. There is also a complex mode with 
qg = 0. The real part is small (Re(e) < 10 x 10~ 13 ) while 
the imaginary part is of order 0(0.01). We ignore them 
as the modes degenerate with the condensate. 



B. Local stability region 

Let us investigate fij for the (1,0,-1), (1,35,0), and 
(1,1,1) vortices. Figures ||(a), ||(b), and ||(c) show the 
local stability regions for each vortex where all the exci- 
tation modes have positive eigenvalues. Figure ^ shows 
the radial shape and qg of the wavefunctions which cor- 
respond to the critical value f2 f . 



1. Similarity to scalar system 

The 0± of the (1,0,-1) vortex for M/N = ±1 and of 
the Alice vortex for M/N = 1 may be understood as 
those of a vortex in a scalar BEC. The in Fig. || re- 
duces to those of the one-component vortex system for 
M/N — 1 because the fully polarized case is nothing but 
a one-component BEC with a single vortex. When we 
decrease M/N towards —1, the (1,0,-1) vortex becomes 
a single-component system again with the opposite wind- 
ing number and the critical O^'s have opposite signs. 

As seen in Fig. |, the modes responsible for f2+ are 
localized on the edge of the condensate (r ~ 3/xm) and 
have a large angular momentum qg ~ 7 when M/N > 0.9 
in both the (1,0,-1) and (1,2,0) vortices. This implies a 
surface-instability mode similar to the surface mode in 
a scalar vortex. The mode remains at the edge of the 
condensate for M/N > 0.8. 

The angular momenta of the modes of Q~ are —1 over 
most of the regions: for < M/N < 1 in the (1,0,-1) 
vortex and for — 1 < M/N < 1 in the Alice vortex. This 
mode is localized at the center when the condensate has 
a large angular momentum [M/N = 1 in Figs. [|(a), ^(b), 
0(c)] . This feature arises from the localized core state in 
a scalar vortex. 



2. Characteristic of spinor system 

The spatial extent of the modes of fl~ , which are lo- 
calized around the core for large M/N expands as M/N 
decreases. This is shown for M/N = in Figs. [I|(a) and 
f|(b), and for M/N = -1 in Figs. ||(c). This is common 
to both the (1,0,-1) and (1,2,0) vortices. The wavefunc- 
tion has a shape similar to that of the condensate and a 
spin structure different from the condensate. Thus this 
instability at fl~ means a spin-flip instability. This is also 
true for n+ of the (1,0,-1) vortex for -0.5 < M/N < 0.5. 
The shape of the wavefunctions is similar to that of the 
condensate. This is not the situation in the antiferro- 
magnetic case. The wavefunction of f2+ indicates a single 
vortex of the -1 component rather than the non- vortex 
-1 component, as shown in Fig. [| (b). 

According to Figs. || (a-c), the critical values f2+ and 
of the (1,0,-1) vortex are close to each other for 
-0.5 < M/N < 0.5. There is no range of ft to stabilize 
the system in the ferromagnetic case because > 
This range is quite narrow in the antiferromagnetic case 
(e.g., 0.261 < < 0.263 for M/N = ±0.2). 

The (1,2,0) vortex reduces to the vortexfree system 
near M/N = -1. As shown in Figs. |(a) and |(b), 
there is an instability for Q — unlike the vortexfree 
scalar system. Figure ^(c) shows the shape of the wave- 
functions u and v of . The shape resembles that of 
the condensate. This instability indicates that the whole 
condensate may have the +1 component without wind- 
ing number. The winding-number combination of the 
(1,2,0) vortex turns into (0,0,0). When // is large and 
negative, 17+ and tt~ become about 0.5 and —0.3 respec- 
tively. Therefore, the system is stable for O = just like 
the vortexfree scalar system. The f2+ shift, depending on 
//, even when M/N has the fixed value —1. This is be- 
cause our calculation does not include the magnetization 
from the non-condensate. 



3. Complex eigenvalues 

Complex eigenvalues emerge when two levels with the 
opposite angular momenta and the same eigenvalues hap- 
pen to appear. The complex eigenvalues of Bogoli- 
ubov equations have been reported for multiply quan- 
tized vortex of scalar BEC |33| and vortex states of bi- 
nary BEC |34| ]. When the eigenvalue e q has a complex 
value, the left-hand side of Eq. ( |l9|) becomes zero. We 
still adopt the real part of such eigenvalues because conti- 
nuity of the real eigenvalues and the real part of complex 
eigenvalues as shown in Fig. [| seems physically meaning- 
ful. 



4. The (1,1,1) vortex 

As for the (1,1,1) vortex, neither the ferromagnetic nor 
the nonmagnetic states can assume intermediate values 
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of magnetization and M/N becomes ±1. Even at each 
magnetization, Slf varies with the chemical potential //. 
The maximum value of Sl~ is about 1.8 (in trap units). 
The wavefunction corresponding to Q,~ is localized at the 
center and has a spin component different from that of 
the condensate. The minimum value of Qj is about 0.4. 
Because Sl~ is always larger, the (1,1,1) vortex is always 
unstable in the sense of local stability. 

This is also the case for the antiferromagnetic system. 
~ 1.8 and £7+ ~ 0.4 mean that the (1,1,1) vortex is 
locally unstable. When the chemical potential pi is large 
enough, these critical values fi^ reduce to 0+ = 0.6 and 
fi~ = 0.3. 



V. DISCUSSION 

We have determined the vortex phase diagram in the 
plane of external rotation frequency SI vs relative mag- 
netization M/N of a spinor BEC system for the cases 
of antiferromagnetic (g s > 0), nonmagnetic (g s = 0), 
and ferromagnetic (g s < 0) interactions. By solving the 
extended Gross-Pitaevskii equation for the spinor F=l 
BEC with the three components (nip = 1,0,-1), we 
have investigated the relative stability of the possible ax- 
isymmetric vortices which are classified according to their 
winding numbers for each component; namely, (1,0,-1), 
(1,2,0), and (1,1,1) for (</>i, 4>q, 4>—i). The excitation spec- 
tra are also studied to find the stability regions for each 
vortex in the SI and M/N plane. This allows us to es- 
timate the vortex nucleation frequency in the actual ex- 
periment (for details, see Refs. j^S], |||] for the correspon- 
dence between local stability and nucleation frequency). 

Yip llill has studied the two non-axisymmetric vortices 
in addition to the present axisymmetric vortices within 
the GP equation under a particular rotation frequency 
(SI ~ 0.4 in our notation) for the antiferromagnetic case 
only. This corresponds to Fig. 2(a) for SI ~ 0.4. Al- 
though the lower M/N region (0 < M/N < 0.2) is 
occupied by a non-axisymmetric vortex, in the remain- 
ing region the (l,x,0) vortex is stabilized over the non- 
axisymmetric vortex, coinciding with our result. We be- 
lieve that the present phase diagram remains valid over a 
broad region even taking into account non-axisymmetric 
vortices. This is because a large part of Yip's phase 
diagram is covered by axisymmetric vortices (see Fig.l 
in Yip's paper fl6f| ). In the future we aim to take into 
account the non-axisymmetric vortex in addition to the 
present axisymmetric ones. This kind of a calculation 
must be performed with the help of the excitation spec- 
trum, which signals the instability towards a more stable 
vortex configuration, for each vortex and is needed to 
predict the experimentally realized vortices in a spinor 
BEC. 




FIG. 4: Spatial variations of ( |u g (r, i) | 2 + \v q (r,i)\ 2 ) as 
a function of the radial direction corresponding to the col- 
lective mode at Sl^ for each M/N is plotted. The (dot- 
ted) lines corresponds to Sl~ (SI*) When the eigenfunctions 
Ui and Vi axe complex number, the sum of the real part 
^2 i (Re(u q (r,i)) 2 + Ke(v q (r, i)) 2 ) is plotted instead. Inset of 
each figure shows the angular momentum index q$ of the mode 
corresponding to SI + . When M/N = 1, the mode of SI t has a 
peak at the edge of the condensate (r ~ 2/im) and the mode at 
Sl~ has a peak at the vortex center (r = 0) . These reproduce 
the result of the scalar BEC. When M/N = 1, qe = 7 and it 
reduces to 1 as M/N decreases, (a) (1,0,-1) vortex in the fer- 
romagnetic case, (b) (1,0,-1) vortex in the antiferromagnetic 
case, (c) (1,1,0) vortex in the nonmagnetic case. 
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